
#include "function.h"
#include <cassert>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
void OutputResult(const VarArray<Variable, 2, 0> &flow,
                  const VarArray<Point, 2, 0> &grid, ostream &os,
                  const vector<double> &positions);
constexpr double re_num{100000};

int main() {
  VarArray<Point, 2, 0> grid;
  GridGenerate(grid, 500, 300);
  OutPutGridInPlot3d(grid);
  // exit(EXIT_FAILURE);
  const int np_in_xdir{grid.Dim(0)};
  const int np_in_ydir{grid.Dim(1)};
  VarArray<Variable, 2, 0> flow(np_in_xdir, np_in_ydir);

  InitializeFlow(flow);

  // 能正常运行， 但 v 不太好， 其它符合的挺好
  CalcUVsecondOrder updatUV(flow, grid, re_num);
  // CalcUVsecondOrder3 updatUV(flow, grid, re_num);

  //   CalcUVsecondOrder2 updatUV(flow, grid, re_num);

  for (int j = 1; j < np_in_xdir; ++j) {
    InitializeNextLine(flow, j);
    // cout << "hello" << endl;
    updatUV.IterateForUV(j);
  }

  ofstream fluid{"flow.plt"};

  vector<double> positions{0.2, 0.6, 1.0};
  OutputResult(flow, grid, fluid, positions);
  fluid.close();

  return 0;
}

void OutputResult(const VarArray<Variable, 2, 0> &flow,
                  const VarArray<Point, 2, 0> &grid, ostream &os,
                  const vector<double> &positions) {
  int npoint_in_ydir = flow.Dim(1);
  int npoint_in_xdir = flow.Dim(0);

  os << "TITLE = \" visualization of the 2Dflow\" \n"
     << "VARIABLES = \"x\", \"y\", \"u\", \"v\"\n";

  os << "ZONE  I=" << npoint_in_xdir << ", J=" << npoint_in_ydir
     << ", DATAPACKING=POINT\n";

  os << std::scientific << std::setprecision(16);

  for (int j = 0; j < npoint_in_ydir; ++j) {
    for (int i = 0; i < npoint_in_xdir; ++i) {
      os << grid(i, j).x << " " << grid(i, j).y << " " << flow(i, j).u << " "
         << flow(i, j).v << "\n";
    }
  }
  // ========================================================================
  // 输出在 1.0 占位处的 blasius 解  positions 中的解
  // string last{"myBlasius.plt"};
  /**
  ofstream blasius{"myBlasius.plt"};
  blasius << std::scientific << std::setprecision(16);
  blasius << "TITLE = \" visualization of the lines\" \n"
          << "VARIABLES = \"my\", \"mu\", \"mv\", \"meta\"\n";

  blasius << "ZONE T=\"myBlasius\",  I=" << npoint_in_ydir << endl;

  for (int j = 0; j < npoint_in_ydir; ++j) {
    const Point &p = grid(npoint_in_xdir - 1, j);
    blasius << p.y << " " << flow(npoint_in_xdir - 1, j).u << " "
            << flow(npoint_in_xdir - 1, j).v << " "
            << p.y / p.x * sqrt(re_num * p.x) << endl;
  }
  blasius.close();
  **/

  ofstream blasius{"myBlasius.plt"};
  blasius << std::scientific << std::setprecision(16);
  blasius << "TITLE = \" visualization of the lines\" \n"
          << "VARIABLES = \"my\", \"mu\", \"mv\", \"meta\"\n";

  for (const auto &pos : positions) {
    string str{"ZONE T=\"myBlasius"};
    // str += to_string(pos);

    double tmp, mmin{1.0e10};
    int index{-1};
    for (int i = 0; i < npoint_in_xdir; ++i) {
      const Point &p = grid(i, 0);
      tmp = fabs(p.x - pos);
      if (tmp < mmin) {
        mmin = tmp;
        index = i;
      }
    }

    str += to_string(grid(index, 0).x);
    blasius << str << "\",  I=" << npoint_in_ydir << endl;

    for (int j = 0; j < npoint_in_ydir; ++j) {
      const Point &p = grid(index, j);
      blasius << p.y << " " << flow(index, j).u << " "
              << flow(index, j).v << " "
              << p.y / p.x * sqrt(re_num * p.x) << endl;
    }
  }

  blasius.close();
  // ==========================================================================
  // 输出边界层厚度  0.99  线性插值
  ofstream thick{"boundary_layer_thickness.plt"};
  thick << std::scientific << std::setprecision(16);
  thick << "TITLE = \" visualization of the boundary layer thickness\" \n"
        << "VARIABLES = \"mx\", \"real_delta\", \"mdelta\"\n";

  thick << "ZONE T=\"mythickness\",  I=" << npoint_in_xdir << endl;

  thick << grid(0, 0).x << "\t" << 0.0 << "\t" << 0.0 << "\n";
  const double delta_u{0.99};
  double y, delt_y, re_x, real_delta;
  for (int i = 1; i < npoint_in_xdir; ++i) {

    re_x = re_num * grid(i, 0).x;
    real_delta = 5.0 * grid(i, 0).x / sqrt(re_x);

    thick << grid(i, 0).x << "\t" << real_delta << "\t";

    y = 0.0;
    for (int j = 0; j < npoint_in_ydir; ++j) {
      if (flow(i, j).u > delta_u) {
        delt_y = grid(i, j).y - grid(i, j - 1).y;
        y = grid(i, j - 1).y + delt_y * ((delta_u - flow(i, j - 1).u) /
                                         (flow(i, j).u - flow(i, j - 1).u));
        thick << y << "\n";
        break;
      }
    }
  }
  thick.flush();
  thick.close();

  // =======================================================================
  ofstream fcf{"theoretical_cf.plt"};
  fcf << std::scientific << std::setprecision(16);
  fcf << "TITLE = \" visualization of the skin friction\" \n"
      << "VARIABLES = \"x\",\"cf\"\n";

  fcf << "ZONE T=\"cf\",  I=" << 4000 << endl;

  double x;
  const double tmpx{1.0e-6};
  double cf;
  for (int i = 1; i <= 1000; ++i) {
    x = i * tmpx;
    cf = 0.664 / sqrt(re_num * x);
    fcf << x << "\t" << cf << "\n";
  }
  const double tmpx2{0.9 / 3000};
  for (int i = 1; i <= 3000; ++i) {
    x = 1.0e-3 + i * tmpx2;
    cf = 0.664 / sqrt(re_num * x);
    fcf << x << "\t" << cf << "\n";
  }
  fcf.flush();
  fcf.close();

  // ======================================================================
  // 二阶近似等距插值
  ofstream fmcf{"my_cf.plt"};
  fmcf << std::scientific << std::setprecision(16);
  fmcf << "TITLE = \" visualization of the skin friction\" \n"
       << "VARIABLES = \"mx\",\"mcf\"\n";

  fmcf << "ZONE T=\"my_cf\",  I=" << npoint_in_xdir - 1 << endl;

  double delty2;
  for (int i = 1; i < npoint_in_xdir; ++i) {
    x = grid(i, 0).x;
    delty2 = grid(i, 2).y - grid(i, 0).y;

    cf = (4.0 * flow(i, 1).u - flow(i, 2).u) / delty2 / re_num * 2.0;

    fmcf << x << "\t" << cf << "\n";
  }
  fmcf.flush();
  fmcf.close();
}